options ls=133;

libname EWORK '.';
libname TU '.';

proc import datafile='new_censusvars_radius.csv' out=newcensus_radius dbms=CSV;
run;

proc import datafile='new_censusvars.csv' out=newcensus dbms=CSV;
run;

proc import datafile='new_censusvars_radius_my.csv' out=newcensus_radius_my dbms=CSV;
run;

proc import datafile='new_censusvars_my.csv' out=newcensus_my dbms=CSV;
run;

*proc means data=newcensus;
*	var family:;
*run;

data newcensus_radius;
	set newcensus_radius;

	lat=floor( lat*1e6 );
	long=floor( long*1e6 );
		
	if pctblack_den>0 then pctblack = pctblack_num/pctblack_den;
	
	if ( income_den>0 ) then do;
		income_lt10 = income_lt10/income_den;
		income_10 = income_10/income_den;
		income_15 = income_15/income_den;
		income_20 = income_20/income_den;
		income_25 = income_25/income_den;
		income_30 = income_30/income_den;
		income_35 = income_35/income_den;
		income_40 = income_40/income_den;
		income_45 = income_45/income_den;
		income_50 = income_50/income_den;
		income_60 = income_60/income_den;
		income_75 = income_75/income_den;
		income_100 = income_100/income_den;
		income_125 = income_125/income_den;
		income_150 = income_150/income_den;
		income_200 = income_200/income_den;
	end;
	else do;
		income_lt10 = .;
		income_10 = .;
		income_15 = .;
		income_20 = .;
		income_25 = .;
		income_30 = .;
		income_35 = .;
		income_40 = .;
		income_45 = .;
		income_50 = .;
		income_60 = .;
		income_75 = .;
		income_100 = .;
		income_125 = .;
		income_150 = .;
		income_200 = .;
	end;
	
	if gt_hs_male_den>0 then gt_hs_male = gt_hs_male_num/gt_hs_male_den;
	if eq_hs_male_den>0 then eq_hs_male = eq_hs_male_num/eq_hs_male_den;
	if gt_hs_female_den>0 then gt_hs_female = gt_hs_female_num/gt_hs_female_den;
	if eq_hs_female_den>0 then eq_hs_female = eq_hs_female_num/eq_hs_female_den;
	
	if publicassistance_den>0 then publicassistance = publicassistance_num/publicassistance_den;
	
	if married_male_den>0 then married_male = married_male_num/married_male_den;
	if married_female_den>0 then married_female = married_female_num/married_female_den;
	if nonmarried_male_den>0 then nonmarried_male = nonmarried_male_num/nonmarried_male_den;
	if nonmarried_female_den>0 then nonmarried_female = nonmarried_female_num/nonmarried_female_den;
	if widowed_male_den>0 then widowed_male = widowed_male_num/widowed_male_den;
	if widowed_female_den>0 then widowed_female = widowed_female_num/widowed_female_den;
	if divorced_male_den>0 then divorced_male = divorced_male_num/divorced_male_den;
	if divorced_female_den>0 then divorced_female = divorced_female_num/divorced_female_den;
	
	if employed_den>0 then employment = employed_num/employed_den;
	if vacant_den>0 then vacant = vacant_num/vacant_den;
	if ownocc_den>0 then ownocc = ownocc_num/ownocc_den;
	if withmort_den>0 then withmort = withmort_num/withmort_den;
	if foreignborn_den>0 then foreignborn = foreignborn_num/foreignborn_den;
	
	if medianrent_a_den>0 then medianrent_a = medianrent_a_num/medianrent_a_den/1000;
	else medianrent_a = .;
	
	if medianrent_b_den>0 then medianrent_b = medianrent_b_num/medianrent_b_den/1000;
	else medianrent_b = .;
	
	if medianhousevalue_den>0 then medianhousevalue = medianhousevalue_num/medianhousevalue_den/1000;
	else medianhousevalue = .;
	
	if mhv_den>0 then mhv = mhv_num/mhv_den/1000;
	else mhv = .;

	if meanincome_den>0 then meanincome = meanincome_num / meanincome_den/1000;
	else meanincome = .;
	
	keep long lat pctblack income_lt10 income_10 income_15 income_20 income_25 income_30 income_35 income_40 income_45 income_50 income_60
	     income_75 income_100 income_125 income_150 income_200 gt_hs_male eq_hs_male gt_hs_female eq_hs_female publicassistance married_male
	     married_female nonmarried_male nonmarried_female widowed_male widowed_female divorced_male divorced_female employment vacant ownocc
	     withmort foreignborn meanincome;
run;

* This should be the same code applied to the modified radius approach (centroid based on census block group internal point);
data newcensus_radius_my;
	set newcensus_radius_my;

	if pctblack_den>0 then pctblack_my = pctblack_num/pctblack_den;
	
	if ( income_den>0 ) then do;
		income_lt10_my = income_lt10/income_den;
		income_10_my = income_10/income_den;
		income_15_my = income_15/income_den;
		income_20_my = income_20/income_den;
		income_25_my = income_25/income_den;
		income_30_my = income_30/income_den;
		income_35_my = income_35/income_den;
		income_40_my = income_40/income_den;
		income_45_my = income_45/income_den;
		income_50_my = income_50/income_den;
		income_60_my = income_60/income_den;
		income_75_my = income_75/income_den;
		income_100_my = income_100/income_den;
		income_125_my = income_125/income_den;
		income_150_my = income_150/income_den;
		income_200_my = income_200/income_den;
	end;
	else do;
		income_lt10_my = .;
		income_10_my = .;
		income_15_my = .;
		income_20_my = .;
		income_25_my = .;
		income_30_my = .;
		income_35_my = .;
		income_40_my = .;
		income_45_my = .;
		income_50_my = .;
		income_60_my = .;
		income_75_my = .;
		income_100_my = .;
		income_125_my = .;
		income_150_my = .;
		income_200_my = .;
	end;
	
	if gt_hs_male_den>0 then gt_hs_male_my = gt_hs_male_num/gt_hs_male_den;
	if eq_hs_male_den>0 then eq_hs_male_my = eq_hs_male_num/eq_hs_male_den;
	if gt_hs_female_den>0 then gt_hs_female_my = gt_hs_female_num/gt_hs_female_den;
	if eq_hs_female_den>0 then eq_hs_female_my = eq_hs_female_num/eq_hs_female_den;
	
	if publicassistance_den>0 then publicassistance_my = publicassistance_num/publicassistance_den;
	
	if married_male_den>0 then married_male_my = married_male_num/married_male_den;
	if married_female_den>0 then married_female_my = married_female_num/married_female_den;
	if nonmarried_male_den>0 then nonmarried_male_my = nonmarried_male_num/nonmarried_male_den;
	if nonmarried_female_den>0 then nonmarried_female_my = nonmarried_female_num/nonmarried_female_den;
	if widowed_male_den>0 then widowed_male_my = widowed_male_num/widowed_male_den;
	if widowed_female_den>0 then widowed_female_my = widowed_female_num/widowed_female_den;
	if divorced_male_den>0 then divorced_male_my = divorced_male_num/divorced_male_den;
	if divorced_female_den>0 then divorced_female_my = divorced_female_num/divorced_female_den;
	
	if employed_den>0 then employment_my = employed_num/employed_den;
	if vacant_den>0 then vacant_my = vacant_num/vacant_den;
	if ownocc_den>0 then ownocc_my = ownocc_num/ownocc_den;
	if withmort_den>0 then withmort_my = withmort_num/withmort_den;
	if foreignborn_den>0 then foreignborn_my = foreignborn_num/foreignborn_den;
	
	*if medianrent_a_den>0 then medianrent_a_my = medianrent_a_num/medianrent_a_den/1000;
	*else medianrent_a_my = .;
	
	*if medianrent_b_den>0 then medianrent_b_my = medianrent_b_num/medianrent_b_den/1000;
	*else medianrent_b_my = .;
	
	*if medianhousevalue_den>0 then medianhousevalue_my = medianhousevalue_num/medianhousevalue_den/1000;
	*else medianhousevalue_my = .;
	
	*if mhv_den>0 then mhv_my = mhv_num/mhv_den/1000;
	*else mhv_my = .;
	
		if meanincome_den>0 then meanincome_my = meanincome_num / meanincome_den/1000;
		else meanincome_my = .;

	keep fips pctblack_my income_lt10_my income_10_my income_15_my income_20_my income_25_my income_30_my income_35_my income_40_my income_45_my income_50_my income_60_my
	     income_75_my income_100_my income_125_my income_150_my income_200_my gt_hs_male_my eq_hs_male_my gt_hs_female_my eq_hs_female_my publicassistance_my married_male_my
	     married_female_my nonmarried_male_my nonmarried_female_my widowed_male_my widowed_female_my divorced_male_my divorced_female_my employment_my vacant_my ownocc_my
	     withmort_my foreignborn_my meanincome_my;
run;

proc sort data=newcensus_radius_my;
	by fips;
run;

data temp;
	set newcensus_my;  * Get population from nearest block group;
	
	if employed_den>0 then population = employed_den;
	else population = .;

	keep fips population;
run;

proc sort data=temp;
	by fips;
run;

* This adds in the population to be used as a weight;
data newcensus_radius_my;
	merge newcensus_radius_my temp;
	by fips;
run;

* Now make a tract-based measure that is a weighted average of the block group results;
data junk;
	set newcensus_radius_my;
	
	fips = floor(fips/10)*10;  * Replace block group with zero;
run;

proc means data=junk noprint;
	by fips;
	weight population;
	where population ne .;
	output out=junk mean=;
run;

data newcensus_radius_my;
	set newcensus_radius_my junk;
run;

data newcensus;
	set newcensus;
	
	lat=floor( lat*1e6 );
	long=floor( long*1e6 );

	if medianrent_a_den>0 then medianrent_a = (medianrent_a_num)/medianrent_a_den/1000;
	else medianrent_a = .;
	
	if medianrent_b_den>0 then medianrent_b = medianrent_b_num/medianrent_b_den/1000;
	else medianrent_b = .;
	
	if medianhousevalue_den>0 then medianhousevalue = medianhousevalue_num/medianhousevalue_den/1000;
	else medianhousevalue = .;

	if mhv_den>0 then mhv = mhv_num/mhv_den/1000;
	else mhv = .;	
	
	if minfips<=0 then minfips = .;
	else minfips = floor( minfips/1e7 );
		
	keep long lat medianrent_a medianrent_b medianhousevalue mhv minfips;
run;

data newcensus_my;
	set newcensus_my;
	
*	lat=floor( lat*1e6 );
*	long=floor( long*1e6 );

	if medianrent_a_den>0 then medianrent_a_my = (medianrent_a_num)/medianrent_a_den/1000;
	else medianrent_a_my = .;
	
	if medianrent_b_den>0 then medianrent_b_my = medianrent_b_num/medianrent_b_den/1000;
	else medianrent_b_my = .;
	
	if medianhousevalue_den>0 then medianhousevalue_my = medianhousevalue_num/medianhousevalue_den/1000;
	else medianhousevalue_my = .;

	if mhv_den>0 then mhv_my = mhv_num/mhv_den/1000;
	else mhv_my = .;	
	
	if minfips<=0 then minfips_my = .;
	else minfips_my = floor( minfips/1e7 );
	
	if employed_den>0 then population = employed_den;  * Should be population 16 or older;
	else population=.;
		
	keep fips medianrent_a_my medianrent_b_my medianhousevalue_my mhv_my minfips_my population;
run;

proc sort data=newcensus_my;
	by fips;
run;

data junk;
	set newcensus_my;
	
	fips = floor(fips/10)*10;
run;

proc means data=junk noprint;
	by fips;
	weight population;
	where population ne .;
	output out=junk sum=;
run;

* This file will have both block group and tract level values;
data newcensus_my;
	set newcensus_my junk;
run;

proc sort data=newcensus_my;
	by fips;
run;

proc sort data=newcensus_radius_my;
	by fips;
run;

data newcensus_my;
	merge newcensus_my newcensus_radius_my;
	by fips;
run;

* Add on the CohenCole replications;
proc sort data=newcensus_radius noduprec;
	by lat long;
run;

proc sort data=newcensus noduprec;
	by lat long;
run;

data newcensus;
	merge newcensus_radius newcensus;
	by lat long;

	counter = 1;
run;

proc sort data=newcensus noduprec;
	by lat long;
run;

data both;
	set EWORK.both;
	where lat ne . and long ne .;
	
	lat = floor( lat*1e6 );
	long = floor(long * 1e6);
run;

proc sort data=both;
	by lat long;
run;

proc sort data=newcensus noduprec;
	by lat long;
run;

data newreplication;
	merge both (in=in1) newcensus (in=in2);
	by lat long;
	
	if in1 & in2;
	
	limit = re28/1000;
	bal = re33/1000;
	util_pct = re34;
	
	if util_pct>=0 then availcredit = limit * ( 1 - util_pct/100 );
	else availcredit = .;
	
	if availcredit ne . then util_dol = limit - availcredit;
	else util_dol = .;
	
	stcnty = state*1000 + county;
	
	if ( (age>0) and (availcredit ne .) and (pctblack ne .) and (medianrent_b>0) and (mhv>0) and (income_lt10>=0) and (tr_am>0) and (stcnty>1000) and (withmort ne .) ) then insample1 = 1;
	else insample1 = 0;
	
	* This one removes the restriction on availcredit;
	if ( (age>0) and (pctblack ne .) and (medianrent_b>0) and (mhv>0) and (income_lt10>=0) and (tr_am>0) and (stcnty>1000) and (withmort ne .) ) then insample2 = 1;
	else insample2 = 0;
	
	* This one removes the restriction on availcredit adds limit restriction;
	if ( (age>0) and (limit>0) and (pctblack ne .) and (medianrent_b>0) and (mhv>0) and (income_lt10>=0) and (tr_am>0) and (stcnty>1000) and (withmort ne .) ) then insample3 = 1;
	else insample3 = 0;
run;

* Now to merge on My method of creating census variables;
proc sort data=newreplication;
	by person year;
run;

data geogdata;
	set TU.kpb_geography2003 (in=in1) TU.kpb_geography2004 (in=in2);
	
	if in1 then year = 2003;
	if in2 then year = 2004;
	
	if ( kpb_state>0 and kpb_county>0 and kpb_tract>0 and kpb_block>0 ) then fips = kpb_state*1e10 + kpb_county*1e7 + kpb_tract*10 + kpb_block;
	else if ( kpb_state>0 and kpb_county>0 and kpb_tract>0 ) then fips = kpb_state*1e10 + kpb_county*1e7 + kpb_tract*10;
	else fips = .;
	
	keep fips person year;
run;

proc sort data=geogdata;
	by fips;
run;

proc sort data=newcensus_my;
	by fips;
run;

data geogdata;
	merge geogdata (in=in1) newcensus_my (in=in2);
	by fips;
	
	if in1 and in2;
run;

proc sort data=geogdata;
	by person year;
run;

data EWORK.newreplication;
	merge newreplication (in=in1) geogdata (in=in2);
	by person year;
	
	stcnty_my = floor(fips/1e7);
	if stcnty_my = 0 then stcnty_my = .;
	
	unusedlim = limit - bal;

	if ( (age>0) and (limit>0) and (pctblack_my ne .) and (medianrent_b_my>0) and (mhv_my>0) and (income_lt10_my>=0) and (tr_am>0) and (stcnty_my>1000) and (withmort_my ne .) ) then insample1_my = 1;
	else insample1_my = 0;
	
	* This one removes the restriction on availcredit;
	if ( (age>0) and (pctblack_my ne .) and (medianrent_b_my>0) and (mhv_my>0) and (income_lt10_my>=0) and (tr_am>0) and (stcnty_my>1000) and (withmort_my ne .) ) then insample2_my = 1;
	else insample2_my = 0;
	
	if in1;
run;

* These procedures generate the data to produce Figure 1;
* Note:  Column (1) is as reported by Cohen-Cole (2011);

* Column (2);
proc means data=EWORK.newreplication n median mean;
	var limit util_pct bal tr_am availcredit util_dol pctblack age
			income_lt10 income_10 income_15 income_20 income_25 income_30 income_35 income_40 income_45 income_50 income_60 income_75 income_100 income_125 income_150 income_200
			foreignborn gt_hs_male gt_hs_female eq_hs_male eq_hs_female publicassistance married_male married_female nonmarried_male nonmarried_female widowed_male widowed_female divorced_male divorced_female
			employment vacant ownocc withmort medianrent_b mhv meanincome;
	where re34 ne .;
	title "RE34 ne .";
run;

* Column (3);
proc means data=EWORK.newreplication n median mean;
	var limit util_pct bal tr_am availcredit util_dol pctblack age
			income_lt10 income_10 income_15 income_20 income_25 income_30 income_35 income_40 income_45 income_50 income_60 income_75 income_100 income_125 income_150 income_200
			foreignborn gt_hs_male gt_hs_female eq_hs_male eq_hs_female publicassistance married_male married_female nonmarried_male nonmarried_female widowed_male widowed_female divorced_male divorced_female
			employment vacant ownocc withmort medianrent_b mhv meanincome;
	where insample1=1;
	title "Insample1=1";
run;

* Column (4);
proc means data=EWORK.newreplication n median mean;
	var limit util_pct bal tr_am availcredit util_dol pctblack age
			income_lt10 income_10 income_15 income_20 income_25 income_30 income_35 income_40 income_45 income_50 income_60 income_75 income_100 income_125 income_150 income_200
			foreignborn gt_hs_male gt_hs_female eq_hs_male eq_hs_female publicassistance married_male married_female nonmarried_male nonmarried_female widowed_male widowed_female divorced_male divorced_female
			employment vacant ownocc withmort medianrent_b mhv meanincome;
	where insample2=1 and insample1=0;
	title "Insample2=1 and insample1=0";
run;

* Not reported;
proc means data=EWORK.newreplication n median mean;
	var limit util_pct bal tr_am availcredit util_dol pctblack age
			income_lt10 income_10 income_15 income_20 income_25 income_30 income_35 income_40 income_45 income_50 income_60 income_75 income_100 income_125 income_150 income_200
			foreignborn gt_hs_male gt_hs_female eq_hs_male eq_hs_female publicassistance married_male married_female nonmarried_male nonmarried_female widowed_male widowed_female divorced_male divorced_female
			employment vacant ownocc withmort medianrent_b mhv meanincome;
	title "All Observations";
run;

endsas;
